%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SVM Decoding
%
% GE_Step1c_average_over_subjects.m
%
% This script will read accuracy results from *.mat files
% for individual subjects, calculate the average and standard
% error of mean (or 95%-confidence intervals), 
% plots them, and saves them into a new *.mat file.
%
% Modified from Dimitrios Pantazis's code 
%  by SA and AS 2018-08-02
%
% Modified from UAG_average_over_subjects 
% by ON 5-AUG-2019
%
% ON & SA 2019-08-09: Modified to averaged over condition pairs 
%                     for each individual subject


%initialize

clear all;
close all;

svm_dir = '/autofs/space/clive_001/users/adriana/GE_SVM';
scripts_dir = sprintf('%s/scripts', svm_dir);
%functions_dir = sprintf('%s/scripts/Functions', svm_dir);

%cd(svm_dir);

addpath(genpath(scripts_dir));
%addpath(genpath(functions_dir));

%SubjectNames = {'Subject01'}; %for more subjects use SubjectNames = {'Subject01' 'Subject02'}; etc

%SubjectNames =     {'GE_05', 'GE_06', 'GE_07', 'GE_08', 'GE_09', 'GE_10', 'GE_11'};
SubjectNames = {'GE_12', 'GE_13', 'GE_14', 'GE_15', 'GE_16'};
n_subjects = length(SubjectNames);

%SubjectDirectories={'GE_05', 'GE_06', 'GE_07', 'GE_08', 'GE_09', 'GE_10', 'GE_11'};
%
%nSubjects = length(SubjectNames);
%
% %Analysis_conditions = {'neighbors_vs_seeds', 'words_vs_seeds', 'nonwords_vs_seeds'};
%Analysis_conditions = {'neighbors_vs_seeds'};
%n_analysis_cond = length(Analysis_conditions);

%analysis_tag = '_neighbors_vs_seeds_16Slices_2Rois_SMG_pMTG'
analysis_tag = ''

% Specify the Regions of Interest ("ROIs", also called "labels":
%roiset_array = {'pMTG'};
roiset_array = {'SMG','pMTG', 'SMG_and_pMTG'};
%roiset_array = {'LIFG'};
n_roisets = length(roiset_array);

% Specify the condition pairs:

train_conditionA_array = {'11','11','11','11','11','22','22','22','22','33','33','33','44','44','55'};
train_conditionB_array = {'22','33','44','55','66','33','44','55','66','44','55','66','55','66','66'};
test_conditionA_array = {  '1', '1', '1', '1', '1', '2', '2', '2', '2', '3', '3', '3', '4', '4', '5'};
test_conditionB_array = {  '2'  '3', '4', '5', '6', '3', '4', '5', '6', '4', '5', '6', '5', '6', '6'};

% train_conditionA_array = {'10','10','10','10','10','20','20','20','20','30','30','30','40','40','50'};
% train_conditionB_array = {'20','30','40','50','60','30','40','50','60','40','50','60','50','60','60'};
% test_conditionA_array  = { '1', '1', '1', '1', '1', '2', '2', '2', '2', '3', '3', '3', '4', '4', '5'};
% test_conditionB_array  = { '2', '3', '4', '5', '6', '3', '4', '5', '6', '4', '5', '6', '5', '6', '6'};

n_condpairs = min(length(train_conditionA_array), length(train_conditionB_array));

% Select the time window for linear fit (in milliseconds):
%fit_tmin = -100;
%fit_tmax = 1000;
fit_tmin = 100;
fit_tmax = 800;

input_fitline='Y';
%input_fitline='N';

%-----------------------------------------------------------------

time_tag = datetime('now', 'Format', 'yyyyMMdd-HHmm');

%for i_analysis_cond = 1:n_analysis_cond
    %analysis_cond = char(Analysis_conditions(i_analysis_cond));
    %resultsdir = sprintf('%s/GE_grand_average_results/%s/%s',svm_dir,analysis_cond,datetime);
    %mkdir(resultsdir);
    
figure('Position',[0, 0, 1650, 1200]); % all subjects and ROIs in one figure

for i_roiset = 1:n_roisets
    roiset = char(roiset_array(i_roiset));

    for i_subject = 1:length(SubjectNames) %for all subjects
        %i_subject = 1;
        subject = char(SubjectNames(i_subject));
        %subject_dir = char(SubjectDirectories(i_subject));
        subject_dir = char(strcat(subject, analysis_tag));

        %inputdir = sprintf('%s/%s_%s/results/', svm_dir, subject, analysis_cond);
        input_dir = sprintf('%s/%s/results/', svm_dir, subject_dir);

        % Create the "<subject>/figures" subdirectory in case it doesn't exist yet
        figure_dir = sprintf('%s/%s/figures', svm_dir, subject_dir);
        [status, msg, msgID] = mkdir(figure_dir);

        %figure(i_subject*100+i_roiset); % separate figure for each ROI (and each subject)
       
        for i_condpair = 1:n_condpairs   
            train_conditionA= char(train_conditionA_array(i_condpair));
            train_conditionB= char(train_conditionB_array(i_condpair));
            test_conditionA= char(test_conditionA_array(i_condpair));
            test_conditionB= char(test_conditionB_array(i_condpair));

            % Get "accuracy" and "Time" from file:
            inputfile = sprintf('%s/%s_%s_train%sand%s_test%svs%s_Accuracy.mat', input_dir, subject, roiset,train_conditionA,train_conditionB,test_conditionA,test_conditionB);
            load(inputfile)

            %plot(Time, accuracy);
            %yline(50);
            %hold on;

            accuracy_allpairs(i_condpair,:) = accuracy; %store decoding accuracy for all condition pairs for this subject
            %accuracy_all(i_subject,:) = accuracy-50;

        end % i_condpair
            
        accuracy_ave = mean(accuracy_allpairs,1);
        accuracy_std = std(accuracy_allpairs,1);
        accuracy_sem = accuracy_std/sqrt(n_condpairs);
        %accuracy_sem = std(accuracy_all,1)/sqrt(nSubjects);
        %conf_level = .95;
        %ts = tinv([1-conf_level conf_level], nSubjects-1);
        %CI = accuracy_ave + ts * accuracy_sem;

        %figure(2);
        %plot(Time, accuracy_ave);
        %hold on;
        %%plot(Time, accuracy_ave-accuracy_std);
        %%plot(Time, accuracy_ave+accuracy_std);
        %yline(50);

        accuracy_allsubjects(i_subject,:) = accuracy_ave;
            
        %resultsfile = sprintf('%s/%s_%s_train%sand%s_test%svs%sAverage_Accuracy.mat',resultsdir, roiset,train_conditionA,train_conditionB,test_conditionA,test_conditionB);
        %save(resultsfile,'accuracy_ave','accuracy_sem','Time','SubjectNames'); % save results

        % Plot the average accuracy:

        subplot(2*n_roisets, n_subjects, (i_roiset-1)*n_subjects +  i_subject);
        plot(Time,accuracy_ave);
        hold on;
        % shading
        plot(Time, accuracy_ave+accuracy_sem);
        plot(Time, accuracy_ave-accuracy_sem);
        %plot(Time, accuracy_ave+accuracy_std);
        %plot(Time, accuracy_ave-accuracy_std);
        hold off;
        yline(50);

        set(gca,'fontsize',8);
        ylabel('Accuracy (%)','fontsize',10)
        xlabel('Time (msec)','fontsize',10)

        %plot_title = subject_dir;
        %maxlength = 33; % cutoff for too long title information 
        %if length(subject_dir) > maxlength
        %    plot_title = subject_dir(1:maxlength);
        %end
        plot_title = subject;
        
        titletext = sprintf('%s %s: all cond. pairs', plot_title, roiset);
        title(titletext, 'fontsize', 10,'Interpreter','none');
        axis([-200 1100 0 100]);
        legend('off');

        %results_fig_file = sprintf('%s/%s/results/%s_%s_%svs%s_AverageAccuracy.png', svm_dir, subject, subject, roiset, condA, condB)
        %print(gcf,[results_fig_file],'-dpng','-r300');
            
    end % i_subject

    accuracy_all_ave = mean(accuracy_allsubjects,1);
    accuracy_all_std = std(accuracy_allsubjects,1);
    accuracy_all_sem = accuracy_all_std/sqrt(n_subjects);
    %accuracy_sem = std(accuracy_all,1)/sqrt(nSubjects);
    %conf_level = .95;
    %ts = tinv([1-conf_level conf_level], nSubjects-1);
    %CI = accuracy_ave + ts * accuracy_sem;

    subplot(2*n_roisets, n_subjects, (n_roisets + i_roiset-1)*n_subjects + 1);
    plot(Time,accuracy_all_ave);
    hold on;
    % shading
    plot(Time, accuracy_all_ave+accuracy_all_sem);
    plot(Time, accuracy_all_ave-accuracy_all_sem);
    %plot(Time, accuracy_all_ave+accuracy_all_std);
    %plot(Time, accuracy_all_ave-accuracy_all_std);
    hold off;
    yline(50);

    set(gca,'fontsize',8);
    ylabel('Accuracy (%)','fontsize',10)
    xlabel('Time (msec)','fontsize',10)

    plot_title = analysis_tag;
    maxlength = 28; % cutoff for too long title information 
    if length(plot_title) > maxlength
        plot_title = plot_title(1:maxlength);
    end
    titletext = sprintf('GE_%s %s: all subjects', plot_title, roiset);
    title(titletext, 'fontsize', 10,'Interpreter','none');
    axis([-200 1100 0 100]);
    legend('off');

end % i_roiset

results_dir = sprintf('%s/GE_grand_average_results/%s/%s',svm_dir,analysis_tag(2:end),time_tag);
[status, msg, msgID] = mkdir(results_dir);

results_fig = sprintf('%s/%s_AverageAccuracy.fig', results_dir,analysis_tag(2:end));
savefig(results_fig);
results_png = sprintf('%s/GE%s_AverageAccuracy.png', results_dir, analysis_tag(2:end));
print(gcf,[results_png],'-dpng','-r300');
    
%end % i_analysis_cond


